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Abstract 

This paper reports a numerical study of the 
Marangoni-Benard (MB) convection in a planar 
fluid layer. The least-squares finite element method 
(LSFEM) is employed to solve the three-dimensional 
Stokes equations and the energy equation. First, 
the governing equations are reduced to be first-order 
by introducing variables such as vorticity and heat 
fluxes. The resultant first-order system is then cast 
into a div-curl-grad formulation, and its elliptic- 
ity and permissible boundary conditions are read- 
ily proved. This numerical approach provides an 
equal-order discretization for velocity, pressure, vor- 
ticity, temperature, and heat conduction fluxes, and 
therefore can provide high fidelity solutions for the 
complex flow physics of the MB convection. Numer- 
ical results reported include the critical Marangoni 
numbers ( M ac ) for the onset of the convection in 
containers with various aspect ratios, and the plan- 
forms of supercritical MB flows. The numerical so- 
lutions compared favorably with the experimental 
results reported by Koschmieder et al.. 
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1 INTRODUCTION 


When a temperature gradient is applied orthogo- 
nally to a thin planar liquid layer with a free inter- 
face, cellular convection occurs from an originally 
quiescent state. The onset of the convection is due 
to the combined effects of the thermal stratification 
instability and the thermo- capillary effect. In par- 
ticular, the temperature dependence of the surface 
tension on the free surface can destabilize the mo- 
tionless fluid state to form regular convective cells. 
Usually, the diameters of these cells are in the same 
order of magnitude as compared to the depth of the 
fluid. This transport phenomenon is referred to as 
the Marangoni-Benard instability due to the first re- 
port of the flow phenomenon by Benard. The name 
also distinguishes it from the Rayleigh-Bernard in- 
stability which could occur without the free surface 
and is induced by the buoyancy. In the past, ex- 
tensive experimental studies of the MB convection 
using silicon oil as the working fluid have been con- 
ducted by Koschmieder et al. [1, 2, 3, 4]. Compre- 
hensive reviews of the MB phenomena can be found 
in Koschmieder [5], Davis [6], and Legros [7]. 

The theoretical studies of the MB convection 
have been focused on the stability analyses. The lin- 
ear stability theory was first established by Pearson 
[8] and later on extended by Nield [9]. Since that 
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time, other type analytical studies flourished, e.g., 
the energy stability theory [10] and the bifurcation 
theory [11, 12]. While these studies have greatly en- 
hanced our understanding of the flow physics, direct 
simulations of the flow phenomena remain attrac- 
tive for further investigation. Full flow equations can 
be numerically solved without assumptions and sim- 
plifications usually employed in the stability anal- 
yses. The complications caused by the buoyancy 
for ground based experiments can also be avoided. 
More than all, the direct simulation is an indispens- 
able tool for studying the regime of the supercritical 
MB flows where few studies have been conducted. 

Duh [13] reported a two-dimensional numeri- 
cal study of the MB flows. A method of stream 
function-vorticity was employed to simulate the MB 
rolls constrained by the bottom and two vertical pe- 
ripheral walls. Numerical results of M ac for the on- 
set of convection as a function of the aspect ratio 
of the container ( A r ) were reported. Here A r is de- 
fined as the ratio of the width of the container to the 
depth. Particularly, he found significant increase of 
Af ac when A r is reduced to be less than two. Win- 
ters et al. [14] used a finite element method to solve 
the two-dimensional flow equations and they inter- 
preted the numerical results by using the bifurcation 
theory. Similar to Duh’s work, they also predicted 
the increase of M ac for lower A r . 

Bestehorn [15] conducted the first three- 
dimensional calculation of the MB convection. A 
special numerical scheme for simulating the three- 
dimensional MB flows was reported. In principle, 
Bestehorn proposed to decompose the divergence 
free velocity into a toroidal and a poloidal parts. For 
fluids with large Prandtl numbers, the toroidal part 
of the velocity is null. As a result, the calculations 
were greatly simplified. By using a spectral method, 
Bestehorn showed the connection between his ampli- 
tude equations and the two-dimensional Ginzburg- 
Landau equation. He presented numerical results of 
the MB planform evolution in containers with very 
large aspect ratios. 

By solving the primitive variables directly, Thess 


et al. [16, 17] reported direct simulations of three- 
dimensional MB flows. Their numerical approach 
took advantage of the flow physics inherent in an 
infinite and periodic MB layer, i.e., the flow motion 
is solely determined by the temperature distribution 
on the free surface. As a result, the calculation pro- 
cedure was simplified and a very efficient pseudo- 
spectral method was developed. The MB flows in 
both weakly and strongly supercritical regimes were 
reported. 

In most practical systems the working fluid is 
bounded by vertical walls, and the wall effects can- 
not be overlooked. For small containers, this situa- 
tion are more pronounced. Rosenblat et al. [18, 19] 
reported the first analytical study of the onset and 
the planform of MB convection in small contain- 
ers. Both linear and non-linear stability analyses 
were conducted. A slippery lateral wall condition 
was employed to avoid the difficulty of the no- slip 
condition. By using a similar analytical method, 
Chen et al. [20] revisited this problem using the 
no- slip condition on the lateral walls. Both studies 
show a sharp increase of M ac as A r decreases below 
2 and no significant increase of M ac for A r > 2. 
Similar conclusion has been reached by Duh us- 
ing two-dimensional direct simulations. Recently, 
Koschmieder and Prahl [3] reported an experimental 
study of the onset and planforms of the MB convec- 
tion in small circular and square containers. This 
study provided the physical evidence to confirm the 
strong increase of the M ac as the A r decreases to a 
small number. In addition, they also reported the 
post-onset Marangoni cells of wedge shapes which 
are not usually seen when using containers of large 
aspect ratios. 

The objective of the present paper is to develop a 
new numerical approach to directly simulate the full 
three-dimensional MB convection. In particular, we 
like to include the no-slip boundary condition on the 
peripheral walls. To this end, we concentrate our at- 
tention to the MB rolls in small square containers. 
As such, we want to recapture the unusual planforms 
observed by Koschmieder and Prahl [3]. Since the 
wall effect must be reckoned, the algorithms used by 
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Bestehorn [15] and Thess et al. [16, 17] cannot be 
employed. To this end, we used the least squares 
finite element method (LSFEM) to solve the equa- 
tions governing the flow physics of the MB convec- 
tion. 

The employed LSFEM is an extension of the work 
originally developed by Jiang et al. [21, 22]. In 
[21], div-curl-grad formulations and their ellipticity 
for incompressible Navier Stokes equations were de- 
rived. In [22], Jiang et al. showed that the LSFEM 
is optimal for the elliptic problems in the sense that 
the global error is of the same order of accuracy as 
compared to the approximation errors. Later on, 
the LSFEM has been extended by Yu et al. [23, 24] 
to solve the compressible viscous flows at low Mach 
numbers. In this work, we shall apply the same tech- 
nique to analyze the ellipticity and the permissible 
boundary conditions for the governing equations of 
the MB convection. 

In Section 2, we present the detailed deriva- 
tion of the first-order formulation for the MB flows, 
including the non-dimensionalization, the div-curl- 
grad system, ellipticity, and the permissible bound- 
ary conditions. In Section 3, the LSFEM and the 
Jacobi conjugate gradient (JCG) method for solv- 
ing the first order equation set are elaborated. In 
the last section, we report the simulated results of 
the MB flows inside square containers. The results 
are compared with the experimental data reported 
by Koschmieder et al. [3]. 

2 THEORETICAL MODELING 

2*1 Governing Equations and Boundary 
Conditions 

In the present study, we like to recapture the the 
MB planforms in square containers reported by 
Koschmieder et al. [3]. In Table 1, which is tab- 
ulated at the end of the text, the properties of the 
silicon oil used in their experiments are listed. Ac- 
cording these data, the Prandtl number of the sil- 


icon oil P T = v/k is about 1000, and the capillary 
number C = pniz/jd is about one thousandth. Note 
that d is the depth of the fluid layer and is in the 
order of mini meter. Since the flow motion is ther- 
mally driven, the Prandtl number is a measure of 
the sluggishness of the fluid velocity; higher Prandtl 
number implies slower motion and vice versa. On 
the other hand, in the absence of gravity the capil- 
lary number C is a measure of the surface deflection. 
And smaller C implies higher surface tension, which 
corresponds to a non-deflecting free surface. Discus- 
sions of the surface deflection effect as a function of 
C can be found in Davis’ work [6]. 

According to the above discussion, two important 
assumptions are made in the present calculation: (1) 
the Prandtl number of the working fluid is large and 
therefore the Stokes equations can be used instead 
of the full Navier Stokes equations, and (2) the cap- 
illary number is small and the free surfaces of the 
MB rolls are flat. As a result, the continuity equa- 
tion, the Stokes equations, and the energy equation 
are considered: 


< 

< 

II 

o 

(i) 

— + v V 2 v = 0, 

9 

(2) 

| + (v-v)r = «v 2 r, 

(3) 


where V = (u, v, w) T is the velocity vector, p is pres- 
sure, p is the density of the fluid, and T is temper- 
ature. The transport properties k and v and the 
density p are assumed constant in the flow field. 

To proceed, the governing equations are reduced 
to a first-order system by introducing new variables: 

n = (t,T),o T = vx v, (4) 

Q = (qx,q y ,qz) T = « V?\ (5) 

where ft is the vorticity with £, r/, and ( as the 
three components, and Q is the heat conduction flux 
vector with q x , q y , and q z as the components in the 
respective directions. This step is necessary for the 
application of the LSFEM so that the C° elements 
can be used in the calculations. As a result, we 
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obtain the following first-order flow equations: 


< 

< 

II 

0 

(6) 

— + v v xfl = 0, 
P 

(7) 

f + CV-V)T=V-Q. 

(8) 


In addition, the vorticity is divergence free and the 
alternative rule of partial differentiation for the heat 
conduction flux vector Q must be satisfied, i.e., 

v * n = o, (9) 

v x Q = 0- (io) 

The boundary conditions on the bottom and the 
side walls of the container is the no-slip condition 
for velocities and vorticities, 

u = v = w = 0 
ft * n = 0 


where n is a unit vector normal to the wall. In addi- 
tion, prescribed temperature at the heated bottom, 
and the insulated condition on the side walls are ap- 
plied. 


T = Tk, on the bottom; 

Q * n = 0, on the vertical walls; 


On the free surface, the Marangoni boundary con- 
ditions are applied. 


du dT 

pv Tz = ~ 7 ^’ 

dv dT 

pv Tz = " 7 V 


(11) 

( 12 ) 


where 7 is the surface tension coefficient. The 
Marangoni conditions represent the relationship be- 
tween the flow shear stress and the tangential sur- 
face tension force across the free surface. Any in- 
homogeneity of the surface tension due to tem- 
perature variations creates a shear force on the 
free surface and therefore results in flow motion. 
These Marangoni boundary conditions are the driv- 
ing force of the Marangoni-Benard convection. By 


using the vorticities and heat conduction fluxes, the 
Marangoni conditions can also be expressed as, 


II 

1 

* I-! 

f 

(13) 

M * 
II 

* 

cl 

(14) 


Since aflat free surface is assumed, we also set w — 0 
on the free surface. 

The heat loss on the free surface is modeled by 
the usual heat transfer condition: 

P KC p ^ = -h(T-T c ) (15) 

where h is a heat transfer coefficient, C v is the con- 
stant pressure specific heat, and T c is the prescribed 
cold temperature of the ambient air. The heat trans- 
fer mechanism on the free surface could be con- 
duction, convection, radiation, and combinations of 
these effects. 


2.2 Non-dimensionalization 

Before the non-dimensionalization, the energy equa- 
tion is reformulated in terms of the temperature per- 
turbation 0, where 


e = T- T ave . ( 16 ) 


A linear distribution of the average temperature T avc 
in the vertical direction is assumed, i.e., 


T ave {z) = T h - Z 1 (T h -T c ), (17) 

where TX and T c are the prescribed hot and cold 
temperatures to set up the MB instability, and d is 
the depth of the liquid layer. This procedure for 
the energy equation is commonly employed in the 
stability analyses of the MB flows. As a result, we 
obtain the following energy equation, 


g + Cv. V )#-=^«v-Q, (i § ) 


where w is the vertical component of the velocity, 
and A T — Th~T c . Note that here we have redefined 
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the heat conduction flux vector Q as the gradient of 
the temperature fluctuation, i.e., Q = k V 9. 


The governing equations and the boundary condi- 
tions are then nondimensionalized by the appropri- 
ate parameters. Here, we choose d and gJ 2 /k as the 
spatial and temporal scales. Therefore, the velocity 
scale is n[d. In addition, the temperature variation 
6 is nondimensionalized by AT. 


u * ud 

c*- 

c ~ V 

p* = 

r pi/K " 

y* — ?- 

~ d \ 

* _ qzd 

Qz - * a : t ’ 


»• = «*, 

- uf. 

*1 K » 

_* JT 

,Y* — 

- *AT’ 
n* 9 

0 ~ AT’ 


^ ~ /c ’ 

?y — «AT’ 


Note that p* is not dimensionless; p* could be inter- 
preted as a dimensionless pressure multiplied by a 
dimensional constant. This treatment is a common 
practice for the Stokes equation. As a result, the 
nondimensionalized flow equations are 


Similarly, after the non-dimensionalization the 
heat convection condition on the free surface be- 
comes, 

q z + Bfi = 0 (29) 

where Bi = hd/npCp is the Biot number, which 
is a dimensionless measure of the heat loss on the 
free surface. Usually, the ambient environment is 
well controlled in the MB experiments and the heat 
transfer on the fluid surface is not efficient. Fur- 
thermore, we note that the energy equation and its 
boundary conditions are formulated in terms of the 
temperature variation 0 instead of temperature it- 
self. Therefore, it is a reasonable assumption to let 
B{ = 0, which implies that the heat transfer to the 
ambient air through the temperature fluctuation 6 
on the free surface is negligible. Instead, all heat 
transfer on the free surface is through the gradient 
of the average temperature T ave . As a result, the 
boundary condition for the energy equation on the 
free surface becomes q z — 0. 


v ■ v = 0 , 

(19) 

VP + VX^ = 0i 

(20) 

| + ( V . V )0-w = vQ, 

(21) 

v ■ fi = o, 

(22) 

v x Q = o, 

(23) 

v x v = n, 

(24) 

= Q 

(25) 


Note that, the superscript * has been dropped for 
convenience. 

The Marangoni boundary conditions are also 
nondimensionalized by the spatial and temporal 
scales, and we get 

tj - - M a q x , (26) 

f = M a q y , (27) 


2.3 Div-Curl-Grad Formulation 

First, we apply the first-order backward differenc- 
ing to the temporal derivative term of the energy 
equation. In addition, by using the definition of the 
heat conduction fluxes, we transform the nonlinear 
convective terms of the energy equation into an alge- 
braic expression. As a result, a new set of first-order 
equations is obtained: 


v-v = o, 

(30) 

VxV=ft, 

(31) 

v • ft = 0 , 

(32) 

V x + VP = °* 

(33) 

v • Q = - s ” _l ) + ( v ■ Q) 

- w, (34) 

v X Q = 0, 

(35) 

V0 = Q> 

(36) 


where 


M a = 


jATd 

p VK 


is the Marangoni number. 


where At is the time step, and the superscript n - 1 
(28) denotes the previous time step. Note that all right 
hand sides are algebraic and they have nothing to 
do with the classification of this equation set. 
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Apparently, these equations are composed of div- 
curl-grad systems: (30) and (31) are a div-curl sys- 
tem for velocities; (32) and (33) are a div-curl system 
for vorticity and pressure; and, (34), (35) and (36) 
are a div-curl-grad system for the heat fluxes and 
temperature fluctuation. As such, we arrive at a 
system with fifteen equations and eleven unknowns. 
The inconsistency between the number of the un- 
knowns and the number of equations has been com- 
monly referred to a s the so-called “over-determined” 
problem. It should be emphasized, however, that 
the “overdetermined” problem is a notion borrowed 
from linear algebra. For partial differential equa- 
tions, this interpretation leads to misconception. In 
the next section, we shall show that the system with 
suitable boundary conditions is elliptic and well- 
posed. 

2.4 Ellipticity 

Since there are eleven unknowns, the equation set 
cannot be classified by the ordinary method, which 
usually requires an even number of unknowns to 
form complex conjugate eigenvalues for elliptic sys- 
tems. To overcome the difficulty, dummy vari- 
ables axe introduced to reconstruct a even-number- 
unknowns system. 

The first-order equation set, (30-36), is actually 
composed of two div-curl systems and a div-curl- 
grad system. In the first system, (30) and (31), the 
div and curl operators operate on three unknowns, 
u, v y and, w , and the system is composed of four 
equations. Here, we introduce a new variable (f>, and 
the system becomes 

V-V = 0, (37) 

V<£+VxV = ft, (38) 

where the dummy variable <p satisfies the boundary 
condition <t> = 0 on T. By applying the divergence 
operator to (38) and considering V ' V X V = 0 and 
y ■ fl = 0 we have 

V 2 <£ = 0 in ft, 

4> = 0 on r. 


Therefore <j> = 0, and the original system has not 
been changed. 

The second div-curl set, (32) and (33), is con- 
structed by four equations and four unknowns i.e., 
f , 17, C, and p, operated by the div and purl opera- 
tors, and no dummy variable is needed. 

The third div-curl-grad set, (34-36), has four un- 
knowns ( 6 , q x , q y , and q z ) determined by seven equa- 
tions. Hence, four dummy variables and one equa- 
tion are introduced into this system: 


v • Q = /» 

(41) 

vtf + V X Q = 0 , 

(42) 

v-* = 0, 

(43) 

yx’f + v^Q, 

(44) 


where t? and = (V’ljtfo,^) are the dummy vari- 
ables, and / is the right hand side of the energy equa- 
tion (see Eq.(34)). As such, this subsystem has eight 
unknowns and eight equations. By applying a diver- 
gence operator to (42) and considering V’VXQ = 0, 
we have V 2 ^ = 0 inside the computational domain. 
Combined with the prescribed boundary condition 
1? = 0 at T, we get t? = 0, and the introduction 
of t? does not change the original equations. Sim- 
ilarly, by applying a curl to (44), and considering 
V X y0 = V x Q = 0, we have 

V x(v x tf) = 0. (45) 

We also know that 

V(VX^) = 0. (46) 

With the boundary condition n x 1 4 r = 0, it can be 
shown that ^ = 0. Therefore, the introduction of 
4 / does not change the original system of equations. 

Facilitated by the dummy variables, we now 
have sixteen equations and sixteen unknowns. In 
the Cartesian coordinates, the first-order system of 
equations have the following form: 

du dv_ dw _ 

Jh + dy + dz ~ ’ 


(39) 

(40) 


6 


(47) 



90 dw 9u _ 

dx + dy dz 

d<j) du dw 

dy dz dx 

90 dv 9u _ 

dz ^ dx dy ’ 

dx + dy T dz ’ 

, ^p = n 
9y dz + 9a: ’ 

_ 0 

<9z <9x <9y 5 

^ _ 0 
dx dy dz ’ 

dqr,dQy,dq 1 _ f 

ax ay oz 

dd_ dq z dqy 

9a: dy dz 

dd dqr dqz 

dy dz dx 

dti dqy dqz 

dz + 9z 9y ’ 
901 902 903 _ 

a I Q T r\ — ”, 

ox oy oz 

d6 903 902 


dx 9y 


9z 


— 


90 90i 903 


9y dz 


dx 


= ?5M 


90 902 90i _ 

9z 9a: 9y ^ 2 ’ 


(48) 

(49) 

(50) 

(51) 

(52) 

(53) 

(54) 

(55) 

(56) 

(57) 

(58) 

(59) 

(60) 
(61) 
(62) 


We then cast the equation set into a matrix form: 


a 4 ? +A 2 ! +a 4 H’ (m) 


in which the unknown vector q is defined as 


q = (u,»,W,0,6^C,P,?r,9y,^,^,01,02,03,7’) T , 

( 64 ) 

and S represents the algebraic terms in the right 
hand side of the equations. The coefficient matrices 
Ai, A 2 , and A3 are tabulated at the end of the text. 


The characteristic polynomial of the system is 
det(AiAi + A2A2 + A3A3) 



c 

0 

0 

0 ^ 

1 

A 

0 

c 

0 

0 


— — G6L 

0 

0 

c 

0 



0 

0 

0 

C ') 


= (detC) 4 






Al 


a 2 

A3 

0 > 

= det 

0 

^3 


-A3 

0 

A2 Aj 
— Ai A2 



> 

Ai 

0 

a 3 j 

— (a? + q + a|) 




^ 0. 

The equation set is indeed elliptic. 


2.5 Permissible Boundary Conditions 


Because the number of equations is even and the 
equation set is elliptic, the required boundary con- 
ditions are standard. In addition, the equation set 
is first-order, therefore only Dirichlet boundary con- 
ditions are permissible. Facilitated by the dummy 
variables, we have sixteen unknowns governed by 
sixteen equations. As such, eight boundary condi- 
tions axe required on each boundary. For the pur- 
pose of discussion, we divide the system of equations 
into two groups: the flow equations, and the heat 
equations. On each boundary, four boundary con- 
ditions are required for each group of the governing 
equations. In Table 2, which is tabulated at the end 
of the text, we propose the permissible boundary 
conditions for the MB convection. 

The outward normal vector of the boundary is 
denoted by n, and the tangential unit vector is 
r = (r!,r 2 ), where r x and r 2 are the two orthog- 
onal components. And Q r i and Q r2 denote the two 
orthogonal components of Q on the free surface. 

For a typical MB flow, there are three type 
boundary conditions: (1) the heated bottom wall 
condition, (2) the insulated vertical wall condition, 
and (3) the free surface condition. For each type 
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boundary, we enforce the known boundary condi- 
tions as listed in Section 2.1. In addition, we in- 
voke the pseudo-boundary conditions of the dummy 
variables to make the proposed system well posed. 
Since the dummy variables are predetermined, the 
adoption of these pseudo-boundary conditions poses 
no theoretical difficulty. These null boundary con- 
ditions are put in parentheses. For the Marangoni 
boundary conditions on the free surface, we specify 
the vorticities by using the heat conduction fluxes 
tangential to the free surface from the previous time 
step. This treatment is in accordance with the flow 
physics that the flow motion is driven by the uneven 
distribution of the surface temperature. Therefore, 
we have proved that for each time step, the flow 
physics described by the semi-discretized governing 
equations is pure elliptic and its solution is deter- 
mined by the abovementioned permissible boundary 
conditions. 


3 THE LEAST-SQUARES FI- 
NITE ELEMENT METHOD 

The LSFEM is used to solve the first order system of 
equations. Due to the flexibility of the LSFEM, the 
number of equations and the number of unknowns 
need not to be identical. Therefore, the dummy vari- 
ables are not included in the numerical solution. The 
first-order system of fifteen equations and eleven un- 
knowns are solved by the LSFEM. A vector form of 
these equations axe considered, 

<«> 

where each entry of the right-hand-side vector S is 
an algebraic equation of eleven dependent variables 
to be solved, i.e., 

Si = Si(qi,q 2 ,--‘,9n) * = 1 , 2,---,15 ( 70 ) 

where q^j = 1, •••,11 are the dependent vari- 
ables. The nonlinear terms are linearized by New- 
ton’s method in the following fashion: 

*r + ' = »r + E (^) cm 


where the superscript m denotes the previous New- 
ton’s method step and m + 1 is the current step. 
Aqj = q™ +l — qj 1 is the increment of the flow vari- 
ables in each Newton’s iteration. After manipula- 
tion, we obtain a new set of equations in vector form 
ready for finite-element discretization, 

» m a , q dq m m 0Aq 

A ° Aq + A ' IT + A ‘ ftT + Ai w + 

A2 dy +Aa dz +A3 dz + j 

To proceed, the governing equations are cast into 
the following operator form: 

LAq = f, (73) 

where the linear operator L is defined as 

L = JC + ^-£ + A?£ + A?£. (74) 

The right-hand-side vector is 

f _ _ A m — ^ - A m — ^ _ A m ^ - S m (75) 
f Al dx Aa dy A * dz { ’ 

We then define the least-squares functional of the 
residual R — LAq — f for admissible Aq as 


J(Aq) = f 
Jn 


R t • R dSl. 


Minimizing the least-squares functional J(Aq) with 
respect to Aq leads to 

6J( Aq) = 0. (77) 


That is, 

[ (LSAqf • (LAq - f) dQ. = 0, (78) 

where 6 denotes the variation of the function. Let 
S Aq = v, and (78) can be written as 

/ (Lv) r (LAq)dft = /(Lvffdfi. (79) 
J n J n 

To employ the finite element method, the compu- 
tational domain is decomposed into N e elements 
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and the element shape functions $,’s are intro- 
duced. The discretized solution in each element 
Aq|(J, J/, z) can be expressed as 

N n 

x, y,z)=Y, $,{x, V, z)(AQ ,-(<))% ( 80 ) 

t=X 

where N n is the number of nodes per element and 
the (AQ t *(f)) c are the nodal values of Aq. The test 
function v is chosen as 

v(x,y,z)=$i(x,y,z) I, ( 81 ) 

where I is the identity matrix. Substituting (80) and 
(81) into (79) gives the linear algebraic equation 

K m AQ = F m , (82) 

where AQ denotes the global nodal values of 

Aq(x,t/), and the final global matrix is 

AT C 

K m = (83) 

e=l 

That is, the global matrix K m is assembled from the 
element matrix (K m ) e , which is defined as 

(Kg)* = / (L*i) T . (»,)■«. (84) 

•'He 

The final right-hand-side vector F m is assembled 
from the element vector (F™) e , which is given as 

(F n e = / (L$i) T -f<W ( 85 ) 

An important feature of the least-squares finite ele- 
ment method can be observed in (79) and (84), i.e., 
the matrix is symmetric. In addition, as long as the 
solution exists, the global matrix is also positive- 
definite. 

The JCG method [25] is employed to invert the 
coefficient matrix. The method is an efficient and 
straight-forward approach for inverting a symmet- 
ric, positive- definite matrix. As long as the solu- 
tion exists, the numerical convergence of the JCG 
method is guaranteed. Because the Jacobi precon- 
ditioning procedure consists of modifying only the 


diagonal terms of the global matrix, the precondi- 
tioned global matrix does not suffer from any fill-in 
and the whole procedure can be implemented in an 
element-by-element fashion such that no global ma- 
trix need to be stored and fine-grain parallelization 
is straightforward. We consider this merit of the LS- 
FEM specially attractive for large scale calculations. 


4 RESULTS AND 

DISCUSSION 

The numerical results reported here are the simu- 
lated MB flows in square containers. The flow fea- 
tures of the MB convection depend on the aspect 
ratio of the containers (A r ) and the Marangoni num- 
ber (Af a ). A r is defined as the ratio of the horizontal 
distance between the opposite walls to the depth of 
the liquid layer. As shown in Fig. 1, the configura- 
tion consists of four insulated vertical walls, a heated 
bottom surface, and a flat free surface. Figure 1 also 
illustrates the specified boundary conditions. For 
the present study, we want to recapture the unusual 
MB planforms reported by Koschmieder et al. [3]. 
Therefore, our attention has been concentrated on 
MB flows in small containers, i.e., A r < 9, where no 
previous numerical study has been reported. Four 
cases are reported: the two, three, four, and five- cell 
MB convection. In all four cases, the same mesh 
(51 X 51 X 19) is used. The mesh is uniform in the 
x and y directions, and is clustered near the free 
surface in the z direction. Although not shown, we 
have conducted the mesh refinement study for the 
four-cell case by doubling the mesh size in each coor- 
dinate direction. Essentially, we have obtained the 
same numerical solution. 

There are three loops of iterations: (1) the outer 
loop is the time marching part; (2) the intermedi- 
ate loop updates the coefficient matrices and source 
terms by Newton’s method; and (3) the inner loop 
solves the variable increment Aq by inverting the 
global matrix using the JCG method. Typically, it 
takes about 50 to 100 time marching steps to con- 
verge about four order of q n — q n_1 , where n de- 
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notes the time step. Further convergence is gener- 
ally much slower. Nevertheless, after about 50 to 
100 time steps, the numerical method usually have 
already caught the MB patterns. Within each time 
step, we perform Newton’s method about three to 
five times to update our coefficient matrix and the 
source vector. For each Newton’s step, it typically 
takes about 300 JCG steps to invert the coefficient 
matrix. 

To start the calculation, we initialize the temper- 
ature field by the initial condition proposed by Thess 
et al. [16]: 

9(x, y, z, 0) = e(x, y)z(2 - z) (86) 

The field e(x, y) is the superposition of all Fourier 
modes supported by the employed numerical mesh. 
The magnitude of e(a:, y ) is set to be one thousandth. 
All other flow properties are initialized by zeros. As 
time evolves, the numerical procedure will pick up 
the most unstable mode and suppress others. In ad- 
dition, we usually start our calculations with very 
low M a and gradually increase the M a until the on- 
set of the flow convection. Finally, the calculation 
will converge to the selected planform. 

Figure 2 shows the numerical solution of a two- 
cell MB convection. Four plots axe shown: (2a) the 
MB planform; (2b) the velocity vectors on the free 
surface; (2c) the temperature contours on the free 
surface; and, (2d) the contours of the vertical vor- 
ticity on the free surface. The planform shown here 
is actually the smoothed contour plot of the veloc- 
ity profile just beneath the free surface. The surface 
topology represents the velocity distribution: the 
bulge-up portion represents the rising flow motion 
and the valley is the downward flow. For M a = 87 
and A t = 5.68, Two triangular MB cells are ob- 
tained. The selected pattern is identical to that re- 
ported by Koschmieder et al. [3]. 

At the centers of the triangles, temperature is 
hotter and therefore the surface tension is lower as 
compared to the area along the walls and the diag- 
onal line, where the temperature is colder. Accord- 
ingly, this unbalanced surface tension force results 


in flow motions from the hot region to the cold re- 
gion; i.e., from the center of the MB cell to the cell 
boundary. To replenish the hot region, hot fluid is 
dragged up from the bottom of the container, and 
therefore the MB convection is sustained. As shown 
in Fig. (2c) The coldest spots on the free surface 
is the higher-left and lower-right corners, where the 
cold fluid is pushed downward to be heated up. It 
is interesting to note that without the effect of the 
vertical walls, the vertical component of the vortic- 
ity (C) is null everywhere. This situation has been 
pointed out by Thess et al. [16]. Figure (2d) shows * 
the distribution of £ on the free surface, and indeed 
all the variations are in the vicinity of the vertical 
walls. This is another indication that our calculation 
has been fairly accurate. 

Figures 3-5 show similar results for three, four, 
and five-cell MB convection for the corresponding 
M a and A T . The patterns are combinations of 
triangles, squares, and wedge-shapes. All these 
patterns have been observed in the laboratory by 
Koschmieder and Prahl [3]. Due to the existence of 
the vertical walls, these patterns are quite different 
from the general conception of the hexagonal MB 
convection usually observed in the containers with 
very large A r . To show the three-dimensional fea- 
tures of the flow fields, we plot vertical sections of 
the temperature profiles and velocity vectors for the 
4-cell case. Figure 6 shows three temperature pro- 
files in the vertical sections: one crosses the center 
of the container, one crosses the center of a MB cell, 
and the last one is close to the wall. For the 4-cell 
case, the coldest spot on the free surface is at the 
center of the container. Figure 7 shows the velocity 
vectors of a vertical section crosses the center of two 
MB cells. It is obvious that there are two upwardly 
rising stream at the centers of the two cells. The 
colder fluid flows downward along the walls and the 
centerline to from several recirculation bubbles. The 
critical Marangoni numbers M ac for all these cases, 
i.e, 5.68 < A r < 8.4, lies between 80 to 85, which 
is consistent with the data reported by Koschmieder 
and Prahl [3]. 
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5 CONCLUDING REMARKS 


pp. 9-15. 


In this paper, we reported the simulations the three- 
dimensional Marangoni-Benard convection based on 
the LSFEM. The continuity equation, the Stokes 
equations, and the time-marching energy equation 
are solved simultaneously. Dependent variables such 
as vorticity and heat conduction fluxes are intro- 
duced to reduce the flow equations to be first-order. 
These first-order equations are composed of sev- 
eral div-curl-grad systems. As such, by using sev- 
eral dummy variables, we show that they are ellip- 
tic. Consequently, the required boundary conditions 
for a well-posed MB flow problem become verita- 
ble. The equation set is solved by the LSFEM, in 
which the coefficient matrix is always symmetric and 
positive- definite. The inversion of the coefficient ma- 
trix is carried out by the JCG method, in which the 
computation is element- by- element and no assem- 
bly of the global matrix is needed. The MB con- 
vections in small square containers with two, three 
four, and five- cell MB convections are simulated. 
The obtained patterns are identical to that reported 
by Koschmieder and Prahl. The critical Marangoni 
numbers for all these cases are also consistent with 
their data. 
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Table 1. Properties of the silicon oil at 25°C 


Symbol 

Property 

Unit 

Value 

V 

Viscosity 

cm 2 /s 

i. 

p 

Density 

OSH 


K 

Thermal diffusivity 

cm 2 /s 

IE 

7 

Surface tension coefficient 

dyne/cm 

13.96 


Table 2 Boundary conditions. 


Conditions 

Flow Eqns. 

Heat Eqns. 

Heated 

Bottom 

u = v = w — 0 
n ■ ft = 0 

9 = 0, Q.i = Q,! = 0 
(9 = 0) 

Insulated 

walls 

o 

II 

if 

Q • n = 0, (n x 'J' = 0) 

Free 

Surface 

n ■ V = 0, Sl r i = M a Q£ L 
fi r2 = -M a C£f l {<f> = 0) 

Q • n = 0, (n x tP = 0) 
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Fig. 1 A schematic of the computational domain. 
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Fig. 2 The Marangoni-Benard convection in a square container with two cells for M a = 87 and 
A r = 5.68. (a) the pattern; (b) the velocity vectors on the free surface; (c) the temperature 
contours on the free surface; and (d) the vertical vorticity contours on the free surface. 




Fig- 3 The Marangoni-Benard convection in a square container with three cells for M a — 95 and 
A t — 6.18. (a) the pattern; (b) the velocity vectors on the free surface; (c) the temperature 
contours on the free surface; and (d) the vertical vorticity contours on the free surface. 







co (d) 


Fig. 4 The Marangoni-Benard convection in a square container with four cells for M a = 95 and 
A r = 6.36. (a) the pattern; (b) the velocity vectors on the free surface; (c) the temperature 
contours on the free surface; and (d) the vertical vorticity contours on the free surface. 



Fig. 5 The Marangoni-Benard convection in a square container with five cells for Af a = 85 and 
A r = 8.48. (a) the pattern; (b) the velocity vectors on the free surface; (c) the temperature 
contours on the free surface; and (d) the vertical vorticity contours on the free surface. 




